MODULE SUVEXCS_MOD
CONTAINS

SUBROUTINE SUVEXCS

USE PARKIND1  , ONLY : JPIM     ,JPRB, JPRD
USE YOMHOOK   , ONLY : LHOOK    ,DR_HOOK, JPHOOK
USE YOS_EXCS  , ONLY : JPRITBL  ,RITBL    ,RIMAX    ,RCHBA    ,&
 & RCHBB    ,RCHBC    ,RCHBD    ,RCHB23A  ,RCHBBCD  ,&
 & RCHBCD   ,RCHETA   ,RCHETB   ,RCHBHDL  ,&
 & RCDHALF  ,RCHBHDL  ,RCDHPI2  ,DRITBL   ,&
 & RLPBB    ,RLPCC    ,RLPDD

! (C) Copyright 1990- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**   *SUBROUTINE* *SUVEXCS* INITIALIZES COMMON BLOCK *YOEVDFS*

!      A. BELJAARS   E.C.M.W.F.   26/03/90

!      PURPOSE
!      -------

!           *SUBROUTINE *SUVEXCS* INITIALIZES THE CONSTANTS NEEDED BY
!      THE EMPIRICAL STABILITY FUNCTIONS AND SETS UP THE TABLE THAT
!      GIVES THE STABILITY PARAMETER ETA (HEIGHT DEVIDED BY
!      *OBUKHOV LENGTH) AS A FUNCTION OF THE GRADIENT *RICHARDSON NUMBER
!      FOR STABLE SITUATIONS.

!      INTERFACE
!      ---------

!           *CALL* *SUVEXCS* FROM *SUSURF*

!      METHOD
!      ------

!           *THE ALGEBRAIC EQUATION TO BE SOLVED IS
!      RI=(PHIH/PHIM**2)*ETA, WHERE *PHIH* AND *PHIM* ARE THE GRADIENT
!      STABILITY FUNCTIONS FOR HEAT AND MOMENTUM RESPECTIVELY.
!           *TO SOLVE THE IMPLICIT ALGEBRAIC EQUATION FOR *ETA* AS A
!      FUNCTION OF *RI*, *NEWTON'S METHOD IS EMPLOYED WITH A FIXED
!      NUMBER OF ITERATIONS. THE NECESSARY DERIVATIVES ARE EVALUATED
!      NUMERICALLY.

!      EXTERNALS
!      ---------

!           *STATEMENT FUNCTIONS *PHIMS*, *PHIHS*

!      Modified

!      P. Viterbo   09/06/2005      P

!      REFERENCE
!      ---------

!           *SEE *PRESS ET AL. (1986; NUMERICAL RECIPES - THE ART OF
!      SCIENTIFIC COMPUTING) FOR DETAILS ON THE SPLINE INTERPOLATION.

!      -----------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM) :: ITMAX, JIT, JJP

REAL(KIND=JPRB) :: ZEPS, ZETA, ZETAM, ZETAP, ZFXM, ZFXP, ZRIB  
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!      INSERT STATEMENT FUNCTIONS (PSI-FUNCTIONS)

#include "fcsvdfs.h"

!*     1. INITIALIZE CONSTANTS IN COMMON BLOCK
!         ---------- --------- -- ------ -----

!     1.1 CONSTANTS RELATED TO ETA(RI)-TABLE

IF (LHOOK) CALL DR_HOOK('SUVEXCS_MOD:SUVEXCS',0,ZHOOK_HANDLE)

RIMAX=10._JPRB
DRITBL  = RIMAX/(JPRITBL-1)
RITBL(1)= 0._JPRB

!     1.2 CONSTANTS FOR THE ELLISON TURNER FUNCTIONS (STABLE)

RCHETA=5._JPRB
RCHETB=4._JPRB
RCHBHDL = 5._JPRB

!     1.21 CONSTANTS FOR HOLTSLAG AND DEBRUIN FUNCTIONS (STABLE PSI)

RCHBA   = 1._JPRB
RCHBB   = 2.0_JPRB/3._JPRB
RCHBC   = 5._JPRB
RCHBD   = 0.35_JPRB
RCHB23A = (2.0_JPRB/3._JPRB)*RCHBA
RCHBBCD = RCHBB*RCHBC/RCHBD
RCHBCD  = RCHBC/RCHBD

!     1.3 CONSTANTS FOR DYER AND HICKS EXPRESSIONS (UNSTABLE)

RCDHALF = 16._JPRB
RCDHPI2 = 2.0_JPRB*ATAN(1.0_JPRB)

!*    3. LOOP OVER TABLE INDEX
!        ---- ---- ----- -----

DO JJP=2, JPRITBL
  ZRIB=(JJP-1)*DRITBL

!     3.1 INITIAL GUESS

  ZETA=RITBL(JJP-1)
  IF (ZRIB  <  0.5_JPRB) THEN
    ITMAX=5
  ELSE
    ITMAX=3
  ENDIF

!     3.2 NEWTON'S ITERATION LOOP WITH DERIVATIVE FROM FINITE DIFFERENCE

  DO JIT=1, ITMAX
    ZEPS=(ZETA+1.0_JPRB)*0.001_JPRB
    ZETAP=ZETA+ZEPS
    ZETAM=ZETA-ZEPS
    ZFXP=ZRIB-ZETAP*PHIHS(ZETAP)/PHIMS(ZETAP)**2
    ZFXM=ZRIB-ZETAM*PHIHS(ZETAM)/PHIMS(ZETAM)**2
    ZETA=ZETA-ZEPS*(ZFXP+ZFXM)/(ZFXP-ZFXM)
  ENDDO

!     3.3 STORE RESULT

  RITBL(JJP)=ZETA

ENDDO

!     CONSTANTS OF THE LOUIS FORMULATION (simplified physics)

RLPBB=5._JPRB
RLPCC=5._JPRB
RLPDD=5._JPRB

IF (LHOOK) CALL DR_HOOK('SUVEXCS_MOD:SUVEXCS',1,ZHOOK_HANDLE)
END SUBROUTINE SUVEXCS
END MODULE SUVEXCS_MOD
